RQ6: original vs (attritioned + imputed)

Description

In this (non-preregistered) analysis, we look at variables collected at years 12 & 16 which have low levels of missingness, and test the effect of selecting only those who took part at year 21 and 26 time points (“attritioning”).

This analysis examines the impact of multiple imputation to correct for attrition bias at three different time points.

We compare ( original ) vs ( attritioned + imputed ) results across 200 multiply imputed datasets.

Load Data

Code
source("0_load_data.R")

df = df %>%
  filter(!(randomfamid %in% exclude_fams_rq1x)) %>%
  filter(!(randomfamid %in% exclude_fams_onesib)) %>%
  filter(!(randomfamid %in% exclude_fams_rq6y))

## range of participation outcomes to use for attrition

range_participation_outcomes = 6:8

# Subset vectors to only include participation outcomes in range
rq1y_twin                         = rq1y_twin[range_participation_outcomes]
rq1y_twin1                        = rq1y_twin1[range_participation_outcomes]
rq1y_twin2                        = rq1y_twin2[range_participation_outcomes]
rq1y_twin_labels                  = rq1y_twin_labels[range_participation_outcomes]
rq1y_twin_labels_clean            = rq1y_twin_labels_clean[range_participation_outcomes]
rq1y_twin_labels_clean_extrashort = rq1y_twin_labels_clean_extrashort[range_participation_outcomes]

boot_compare_results = readRDS(file.path("results","8_rq6_boot_compare_results.Rds"))

# Number of imputations per timepoint
n_imputations = length(boot_compare_results)/length(range_participation_outcomes)

Focal Variables and timepoints

Attrition time points

Code
cbind(
  rq6y,
  rq6y_labels
) %>%
  knitr::kable(caption = "Outcome Variables")
Outcome Variables
rq6y rq6y_labels
lcmfqt1 Y12: Depression (MFQ)
lsdqext1 Y12: Externalising
lcg1 Y12: Cognitive ability
pcexgcsecoregrdm1 Y16: GCSE core subjects grade
Code
cbind(
  rq1y_twin1,
  rq1y_twin_labels_clean
) %>%
  knitr::kable(caption = "Attrition time points")
Attrition time points
rq1y_twin1 rq1y_twin_labels_clean
u1cdata1 Y21 (TEDS21 phase-1 questionnaire)
zmhdata1 Y26 (TEDS26 questionnaire)
zcdata1 Y26 (CATSLife web tests)

Table of all variables used in multiple imputation

Code
# Function to map variable prefix to study wave
get_study_wave = function(var_name) {
  prefix = substr(var_name, 1, 1)
  wave_map = c(
    "a" = "1st Contact",
    "b" = "2 Year",
    "c" = "3 Year",
    "d" = "4 Year",
    "e" = "In Home",
    "g" = "7 Year",
    "h" = "8 Year",
    "i" = "9 Year",
    "j" = "10 Year",
    "l" = "12 Year",
    "n" = "14 Year",
    "p" = "16 Year",
    "r" = "18 Year",
    "u" = "21 Year",
    "z" = "26 Year"
  )
  return(wave_map[prefix])
}

rq6_imputation_vars = c(rq6y_all, rq6z_vars)
impute_vars_labels = var_to_label(rq6_imputation_vars) %>%
  sapply(., function(x) ifelse(is.null(x[1]), "", x[1]))

# Create formatted table for imputation variables
v_impute = data.frame(
  Description = impute_vars_labels,
  `Teds Code` = ifelse(rq6_imputation_vars %in% original_colnames, rq6_imputation_vars, paste0(rq6_imputation_vars, "*")),
  `Range or Level` = sapply(rq6_imputation_vars, function(var) {
    if (var %in% colnames(df)) {
      if (class(df[[var]]) == "numeric") {
        paste0(round(min(df[[var]], na.rm = TRUE), 2), " - ", round(max(df[[var]], na.rm = TRUE), 2))
      } else if (is.factor(df[[var]])) {
        factor_levels = levels(df[[var]])
        paste(c(paste0(factor_levels[1], "*"), factor_levels[-1]), collapse = ", ")
      } else {
        paste(unique(df[[var]]), collapse = ", ")
      }
    } else {
      "Variable not found"
    }
  }),
  N = sapply(rq6_imputation_vars, function(var) {
    if (var %in% colnames(df)) {
      sum(!is.na(df[[var]]))
    } else {
      0
    }
  }),
  `Study.Wave` = sapply(rq6_imputation_vars, get_study_wave)
)

gt(v_impute) %>%
  cols_label_with(fn = ~ gsub("\\.", " ", .x)) %>%
  tab_style(
    style = cell_text(size = px(12)),
    locations = cells_body()
  ) %>%
  tab_style(
    style = cell_text(size = px(12)),
    locations = cells_column_labels()
  ) %>%
  cols_hide(N) %>%
  tab_source_note(
    source_note = "Note: Range or Level shows min-max values for numeric variables or factor levels for categorical variables (reference level marked with *). Variable codes with an asterisk (*) have been derived or modified from the original dataset."
  ) %>%
  tab_options(
    table.width = "70%"
  ) %>%
  cols_width(
    Description ~ px(300),
    everything() ~ px(120)
  )
Description Teds Code Range or Level Study Wave
MFQ scale from 11 MFQ items (child self-report) at 12, 0-22 lcmfqt1 0 - 22 12 Year
SDQ Externalising scale at 12 lsdqext1* 0 - 20 12 Year
G composite scale from child web tests at 12, standardised lcg1 -3.67 - 3.04 12 Year
Core subjects (English, maths, science): mean grade in GCSE results (twin exams at 16), 4-11 pcexgcsecoregrdm1 4 - 11 16 Year
MFQ scale from 11 MFQ items (child self-report) at 12, 0-22 lcmfqt2 0 - 22 12 Year
lsdqext2* 0 - 20 12 Year
G composite scale from child web tests at 12, standardised lcg2 -3.67 - 3.04 12 Year
Core subjects (English, maths, science): mean grade in GCSE results (twin exams at 16), 4-11 pcexgcsecoregrdm2 4 - 11 16 Year
School cohort, see value labels cohort Cohort 1: twins born Jan-94 to Aug-94*, Cohort 2: twins born Sep-94 to Aug-95, Cohort 3: twins born Sep-95 to Aug-96, Cohort 4: twins born Sep-96 to Dec-96 3 Year
Age in years of natural mother at time of birth of twins amumagetw 16 - 45 1st Contact
Single Parent asingle* cohabiting biological mother and father / cohabiting biological parent with other*, single parent 1st Contact
Mother SOC employment level (1st Contact), 1-9 amosoc2* caring for children at home*, 1, 2, 3, 4, 5, 6, 7, 8, 9, no job 1st Contact
Maternal Education (formatted as numeric) amohqualn* 1 - 8 1st Contact
Ethnic origin of twins, simplified coding (1st Contact), 1=white 0=other aethnic 0 - 1 1st Contact
Main language spoken at home (1st Contact), see value labels alang other*, English, English + other 1st Contact
Number of older siblings (formatted as numeric variable) anoldsibn* 0 - 5 1st Contact
Number of younger siblings (formatted as numeric variable) anyngsibn* 0 - 2 1st Contact
Member of a Twins Club (1st Contact), 1Y 0N atwclub 0*, 1 1st Contact
Census data 2001 (code KS001) linked to 1998 postcode: population density, N per hectare cens01pop98density 0 - 310 3 Year
Smoked cigarettes while pregnant (1st Contact), 1Y 0N asmoke 0*, 1 1st Contact
G: overall cognitive ability composite (4 Year), standardised drawg1 -4.47 - 2.62 4 Year
G: overall cognitive ability composite (4 Year), standardised drawg2 -4.47 - 2.62 4 Year
SDQ Conduct scale (4 Year), 0-10 dsdqcont1 0 - 10 4 Year
SDQ Conduct scale (4 Year), 0-10 dsdqcont2 0 - 10 4 Year
SDQ Hyperactivity scale (4 Year), 0-10 dsdqhypt1 0 - 10 4 Year
SDQ Hyperactivity scale (4 Year), 0-10 dsdqhypt2 0 - 10 4 Year
ARBQ total scale (4 Year), 0-36 danxt1 0 - 31 4 Year
ARBQ total scale (4 Year), 0-36 danxt2 0 - 31 4 Year
Chaos overall scale (4 Year), standardised dchatot -2.24 - 4.04 4 Year
Child G composite (7 year twin phone), standardised gcg1 -4.36 - 4.84 7 Year
Child G composite (7 year twin phone), standardised gcg2 -4.36 - 4.84 7 Year
2-subject (maths, English) mean NC level (7 year teacher), 0-4 gt2ac1 0 - 4 7 Year
2-subject (maths, English) mean NC level (7 year teacher), 0-4 gt2ac2 0 - 4 7 Year
SDQ hyperactivity total (7 year parent), 0-10 gpsdqhypt1 0 - 10 7 Year
SDQ hyperactivity total (7 year parent), 0-10 gpsdqhypt2 0 - 10 7 Year
SDQ conduct total (7 year parent), 0-10 gpsdqcont1 0 - 10 7 Year
SDQ conduct total (7 year parent), 0-10 gpsdqcont2 0 - 10 7 Year
SDQ emotion total (7 year parent), 0-10 gpsdqemot1 0 - 10 7 Year
SDQ emotion total (7 year parent), 0-10 gpsdqemot2 0 - 10 7 Year
ARBQ overall anxiety total (7 year parent), 0-52 gpanxt1 0 - 46 7 Year
ARBQ overall anxiety total (7 year parent), 0-52 gpanxt2 0 - 46 7 Year
Conners ADHD overall Total at 8 (0-54) hconnt1 0 - 54 8 Year
Conners ADHD overall Total at 8 (0-54) hconnt2 0 - 54 8 Year
G composite (9 year child), standardised icg1 -3.63 - 2.24 9 Year
G composite (9 year child), standardised icg2 -3.63 - 2.24 9 Year
3-subject (English, maths, science) mean NC level (9 year teacher), 1-5 it3ac1 1 - 5 9 Year
SDQ Conduct scale (child self-report) at 9, 0-10 icsdqcont1 0 - 10 9 Year
SDQ Conduct scale (child self-report) at 9, 0-10 icsdqcont2 0 - 10 9 Year
SDQ Emotion scale (child self-report) at 9, 0-10 icsdqemot1 0 - 10 9 Year
SDQ Emotion scale (child self-report) at 9, 0-10 icsdqemot2 0 - 10 9 Year
SDQ Hyperactivity scale (child self-report) at 9, 0-10 icsdqhypt1 0 - 10 9 Year
SDQ Hyperactivity scale (child self-report) at 9, 0-10 icsdqhypt2 0 - 10 9 Year
SDQ Hyperactivity scale (parent) at 9, 0-10 ipsdqhypt1 0 - 10 9 Year
SDQ Hyperactivity scale (parent) at 9, 0-10 ipsdqhypt2 0 - 10 9 Year
SDQ Conduct scale (parent) at 9, 0-10 ipsdqcont1 0 - 10 9 Year
SDQ Conduct scale (parent) at 9, 0-10 ipsdqcont2 0 - 10 9 Year
SDQ Emotion scale (parent) at 9, 0-10 ipsdqemot1 0 - 10 9 Year
SDQ Emotion scale (parent) at 9, 0-10 ipsdqemot2 0 - 10 9 Year
Negative parental feelings scale, child-self-rated at 9, 0-8 icparnegt1 0 - 8 9 Year
Negative parental feelings scale, child-self-rated at 9, 0-8 icparnegt2 0 - 8 9 Year
ARBQ overall anxiety total (9 year parent), 0-50 ipanxt1 0 - 38 9 Year
ARBQ overall anxiety total (9 year parent), 0-50 ipanxt2 0 - 38 9 Year
Parent Chaos scale at 9, 0-12 ipchatot 0 - 12 9 Year
Child G composite (10 year twin web), standardised jcg1 -3.67 - 2.98 10 Year
Child G composite (10 year twin web), standardised jcg2 -3.67 - 2.98 10 Year
3-subject (English, maths, science) mean NC level (10 year teacher), 1-5 jt3ac1 1 - 5 10 Year
3-subject (English, maths, science) mean NC level (10 year teacher), 1-5 jt3ac2 1 - 5 10 Year
G composite scale from child web tests at 14, standardised ncg1 -4.12 - 3.17 14 Year
G composite scale from child web tests at 14, standardised ncg2 -4.12 - 3.17 14 Year
End of KS3 3-subject Academic achievement mean level (from parent SLQ), 1-9 npks3t3a1 1 - 9 14 Year
End of KS3 3-subject Academic achievement mean level (from parent SLQ), 1-9 npks3t3a2 1 - 9 14 Year
Conners Hyperactivity-Impulsivity scale at 14 (child), 0-27 ncconhit1 0 - 27 14 Year
Conners Hyperactivity-Impulsivity scale at 14 (child), 0-27 ncconhit2 0 - 27 14 Year
Conners total scale at 14 (parent), 0-54 npconnt1 0 - 54 14 Year
Conners total scale at 14 (parent), 0-54 npconnt2 0 - 54 14 Year
Negative Parental Feelings scale at 14 (child), 0-8 ncparnegt1 0 - 8 14 Year
Negative Parental Feelings scale at 14 (child), 0-8 ncparnegt2 0 - 8 14 Year
Chaos at home total scale at 14 (child), 0-12 ncchato1 0 - 11 14 Year
Chaos at home total scale at 14 (child), 0-12 ncchato2 0 - 11 14 Year
G composite scale from child web tests at 16, standardised pcg1 -2.86 - 4.06 16 Year
G composite scale from child web tests at 16, standardised pcg2 -2.86 - 4.06 16 Year
MFQ total scale (child behaviour qnr at 16), 0-26 pcbhmfqt1 0 - 26 16 Year
MFQ total scale (child behaviour qnr at 16), 0-26 pcbhmfqt2 0 - 26 16 Year
SDQ Emotion scale (child behaviour qnr at 16), 0-10 pcbhsdqemot1 0 - 10 16 Year
SDQ Emotion scale (child behaviour qnr at 16), 0-10 pcbhsdqemot2 0 - 10 16 Year
SDQ Conduct scale (child behaviour qnr at 16), 0-10 pcbhsdqcont1 0 - 10 16 Year
SDQ Conduct scale (child behaviour qnr at 16), 0-10 pcbhsdqcont2 0 - 10 16 Year
SDQ Hyperactivity scale (child behaviour qnr at 16), 0-10 pcbhsdqhypt1 0 - 10 16 Year
SDQ Hyperactivity scale (child behaviour qnr at 16), 0-10 pcbhsdqhypt2 0 - 10 16 Year
SDQ Conduct scale (parent behaviour qnr at 16), 0-10 ppbhsdqcont1 0 - 10 16 Year
SDQ Conduct scale (parent behaviour qnr at 16), 0-10 ppbhsdqcont2 0 - 10 16 Year
SDQ Hyperactivity scale (parent behaviour qnr at 16), 0-10 ppbhsdqhypt1 0 - 10 16 Year
SDQ Hyperactivity scale (parent behaviour qnr at 16), 0-10 ppbhsdqhypt2 0 - 10 16 Year
Chaos total score (child web at 16), 0-12 pcchatot1 0 - 12 16 Year
Chaos total score (child web at 16), 0-12 pcchatot2 0 - 12 16 Year
UCAS point score total estimate from known UK exam results (twin qnr at 18) rcqucast1 20 - 910 18 Year
UCAS point score total estimate from known UK exam results (twin qnr at 18) rcqucast2 20 - 910 18 Year
Twin probable highest level of qualification including current study (TEDS21 phase 1 twin qnr), 1-11 see value labels u1chqualp1 1 - 11 21 Year
Twin probable highest level of qualification including current study (TEDS21 phase 1 twin qnr), 1-11 see value labels u1chqualp2 1 - 11 21 Year
SDQ Conduct total score (TEDS21 phase 1 twin qnr), 0-10 u1csdqcont1 0 - 9 21 Year
SDQ Conduct total score (TEDS21 phase 1 twin qnr), 0-10 u1csdqcont2 0 - 9 21 Year
SDQ Hyperactivity total score (TEDS21 phase 1 twin qnr), 0-10 u1csdqhypt1 0 - 10 21 Year
SDQ Hyperactivity total score (TEDS21 phase 1 twin qnr), 0-10 u1csdqhypt2 0 - 10 21 Year
Demographics item: highest qualification ordinal level (TEDS26 twin MHQ), see value labels zmhhqual1 1 - 11 26 Year
Demographics item: highest qualification ordinal level (TEDS26 twin MHQ), see value labels zmhhqual2 1 - 11 26 Year
SDQ Conduct total score (TEDS26 twin MHQ), 0-10 zmhsdqcont1 0 - 9 26 Year
SDQ Conduct total score (TEDS26 twin MHQ), 0-10 zmhsdqcont2 0 - 9 26 Year
SDQ Hyperactivity total score (TEDS26 twin MHQ), 0-10 zmhsdqhypt1 0 - 10 26 Year
SDQ Hyperactivity total score (TEDS26 twin MHQ), 0-10 zmhsdqhypt2 0 - 10 26 Year
Note: Range or Level shows min-max values for numeric variables or factor levels for categorical variables (reference level marked with *). Variable codes with an asterisk (*) have been derived or modified from the original dataset.

Data Processing

Transform Raw Results to Dataframe

Code
# Create pairwise correlation variable names (if needed for cor_resid later)
test_correlation_matrix = matrix(
  nrow = length(rq6y),
  ncol = length(rq6y)
)

for(i in seq_along(rq6y)){
  for(j in seq_along(rq6y)){
    test_correlation_matrix[i,j] = paste(rq6y[i], rq6y[j], collapse = " ")
  }
}

x = test_correlation_matrix[lower.tri(test_correlation_matrix, diag = FALSE)]

x_var = str_extract(x, "^\\S+")
y_var = str_extract(x, "\\S+$")

rm(test_correlation_matrix,x)

# Split results by timepoint
timepoint_names = rq1y_twin
imputed_comparisons_df = list()

for(tp in seq_along(timepoint_names)){

  timepoint_name = timepoint_names[tp]

  # Get all imputations for this timepoint
  tp_pattern = paste0("^", timepoint_name)
  tp_indices = grep(tp_pattern, names(boot_compare_results))

  # Extract results for each imputation
  tp_results = boot_compare_results[tp_indices]

  # For each metric (md, smd, var), aggregate across imputations
  # Using Rubin's rules for multiple imputation

  metrics = c("md", "smd", "var")
  metric_results = list()

  for(metric_idx in seq_along(metrics)){

    metric = metrics[metric_idx]
    
    metric_data = lapply(tp_results, function(imp) {
      # imp[[metric]] is a list of bootstrap iterations
      # Convert to matrix where rows = bootstrap iterations, cols = variables
      sapply(imp[[metric]], function(boot_iter) boot_iter)
    })

    all_boots = do.call(cbind, metric_data)

    # Calculate summary statistics across all bootstrap samples
    metric_summary = apply(all_boots, 1, function(xx) .mean_qi_pd(xx))
    metric_summary_df = do.call(rbind, metric_summary)

    metric_summary_df = data.frame(metric_summary_df)
    metric_summary_df$dataset = rq1y_twin_labels_clean_extrashort[tp]
    metric_summary_df$stat = metric
    metric_summary_df$variable = rq6y_labels
    rownames(metric_summary_df) = NULL

    metric_results[[metric_idx]] = metric_summary_df
  }

  # Combine metrics for this timepoint
  imputed_comparisons_df[[tp]] = do.call(rbind.data.frame, metric_results)
}

# Combine all timepoints into single dataframe
imputed_comparisons_df0 = do.call(rbind.data.frame, imputed_comparisons_df)
rownames(imputed_comparisons_df0) = NULL

# Save for future use
saveRDS(imputed_comparisons_df0, file.path("results", "8_imputed_comparisons_df.Rds"))

Create Table

Code
# Process the dataframe
# Note: Variance values are already percentage change from compare_var function
imputed_comparisons_df = imputed_comparisons_df0 %>%
  mutate(
    variable = rq6y_labels[match(.$variable, rq6y_labels)]
  ) %>%
  group_by(dataset, stat) %>%
  mutate(
    pval_adj = stats::p.adjust(pval, method = "holm")
  )

Means, SMDs, Variances

GT Results Table

Code
imputed_comparisons_df %>%
  filter(stat %in% c("md", "smd", "var")) %>%
  select(-starts_with("."),-pd) %>%
  pivot_wider(
    names_from = "stat",
    values_from = c("y","ymin","ymax","pval","pval_adj"),
    id_cols = c("variable", "dataset")
  ) %>%
  dplyr::rename(Outcome = "variable") %>%
  gt(groupname_col = "dataset") %>%
  tab_options(
    table.width = px(800)
  ) %>%

  # Create column spanners (nested headers)
  tab_spanner(
    label = "Mean Difference",
    columns = c(y_md, ymin_md, ymax_md, pval_md, pval_adj_md)
  ) %>%
  tab_spanner(
    label = "Standardized Mean Difference",
    columns = c(y_smd, ymin_smd, ymax_smd, pval_smd, pval_adj_smd)
  ) %>%
  tab_spanner(
    label = "Variance % Change",
    columns = c(y_var, ymin_var, ymax_var, pval_var, pval_adj_var)
  ) %>%

  # Rename columns under each spanner
  cols_label(
    y_md = "Est", ymin_md = "LB", ymax_md = "UB", pval_md = md("p<sub>unadj</sub>"), pval_adj_md = md("p<sub>adj</sub>"),
    y_smd = "Est", ymin_smd = "LB", ymax_smd = "UB", pval_smd = md("p<sub>unadj</sub>"), pval_adj_smd = md("p<sub>adj</sub>"),
    y_var = "Est", ymin_var = "LB", ymax_var = "UB", pval_var = md("p<sub>unadj</sub>"), pval_adj_var = md("p<sub>adj</sub>")
  ) %>%

  # Format numbers
  fmt(
    columns = !contains("var") & !contains("Outcome"),
    fns = function(x) {gbtoolbox::apa_num(x, n_decimal_places = 3)}
  ) %>%
  fmt(
    columns = c("pval_var", "pval_adj_var"),
    fns = function(x) {gbtoolbox::apa_num(x, n_decimal_places = 3)}
  ) %>%
  fmt_percent(
    columns = c("y_var", "ymin_var","ymax_var"),
    decimals = 2,
    drop_trailing_zeros = FALSE,
    drop_trailing_dec_mark = FALSE
  ) %>%
  # Styling - uniform font size
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_column_spanners()
  ) %>%
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_body()
  ) %>%
  tab_style(
    style = cell_text(size = px(10), align = "center"),
    locations = cells_row_groups()
  ) %>%
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_column_labels()
  ) %>%
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_footnotes()
  ) %>%

  # Highlight significant results - positive effects (light green)
  tab_style(
    style = cell_fill(color = "#d5e8d4"),
    locations = cells_body(
      columns = c(y_md, ymin_md, ymax_md, pval_md, pval_adj_md),
      rows = pval_adj_md < 0.05 & y_md > 0
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#d5e8d4"),
    locations = cells_body(
      columns = c(y_smd, ymin_smd, ymax_smd, pval_smd, pval_adj_smd),
      rows = pval_adj_smd < 0.05 & y_smd > 0
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#d5e8d4"),
    locations = cells_body(
      columns = c(y_var, ymin_var, ymax_var, pval_var, pval_adj_var),
      rows = pval_adj_var < 0.05 & y_var > 0
    )
  ) %>%

  # Highlight significant results - negative effects (light red)
  tab_style(
    style = cell_fill(color = "#f8cecc"),
    locations = cells_body(
      columns = c(y_md, ymin_md, ymax_md, pval_md, pval_adj_md),
      rows = pval_adj_md < 0.05 & y_md < 0
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#f8cecc"),
    locations = cells_body(
      columns = c(y_smd, ymin_smd, ymax_smd, pval_smd, pval_adj_smd),
      rows = pval_adj_smd < 0.05 & y_smd < 0
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#f8cecc"),
    locations = cells_body(
      columns = c(y_var, ymin_var, ymax_var, pval_var, pval_adj_var),
      rows = pval_adj_var < 0.05 & y_var < 0
    )
  ) %>%

  # Add footnotes
  tab_footnote(
    footnote = md("<em>Note.</em> Est = Estimate, LB = Lower Bound 95% Confidence Interval, UB = Upper Bound 95% Confidence Interval. Significant (p<sub>Bonferroni-Holm</sub>) effects are highlighted in green (increases) or red (decreases)."),
    placement = "right"
  ) %>%
  tab_footnote(
    footnote = "P values are Bonferroni-Holm adjusted within each attrition timepoint",
    locations = cells_column_labels(columns = contains("pval_adj")),
    placement = "right"
  ) %>%
  tab_footnote(
    footnote = md("Standardized mean difference = (Mean<sub>attrition + imputed</sub> - Mean<sub>original</sub>) / SD<sub>original</sub>"),
    locations = cells_column_spanners(spanners = "Standardized Mean Difference"),
    placement = "right"
  ) %>%
  tab_footnote(
    footnote = md("Percentage change = (Var<sub>attrition + imputed</sub> - Var<sub>original</sub>) / Var<sub>original</sub>"),
    locations = cells_column_spanners(spanners = "Variance % Change"),
    placement = "right"
  ) %>%
  opt_footnote_marks(marks = c("*", "†","‡"))
Outcome
Mean Difference
Standardized Mean Difference*
Variance % Change†
Est LB UB punadj padj‡ Est LB UB punadj padj‡ Est LB UB punadj padj‡
Y21
Y12: Depression (MFQ) .141 .086 .194 .000 .000 .042 .026 .058 .000 .000 10.05% 5.38% 15.20% .000 .000
Y12: Externalising .036 -.015 .087 .172 .220 .010 -.004 .025 .172 .220 2.85% 0.17% 5.70% .034 .102
Y12: Cognitive ability .011 -.003 .024 .110 .220 .011 -.003 .024 .110 .220 1.04% −1.30% 3.44% .395 .790
Y16: GCSE core subjects grade .012 .001 .024 .033 .100 .010 .001 .019 .033 .100 −0.12% −1.85% 1.69% .865 .865
Y26 (q'aire)
Y12: Depression (MFQ) .172 .110 .231 .000 .000 .051 .033 .069 .000 .000 10.77% 5.40% 16.38% .000 .000
Y12: Externalising .049 -.011 .103 .097 .146 .014 -.003 .030 .097 .146 5.08% 2.00% 8.43% .001 .004
Y12: Cognitive ability .015 -.001 .032 .073 .146 .015 -.001 .032 .073 .146 1.73% −0.88% 4.44% .198 .198
Y16: GCSE core subjects grade .022 .009 .034 .000 .001 .018 .007 .028 .000 .001 −1.76% −3.64% 0.10% .066 .132
Y26 (web test)
Y12: Depression (MFQ) .383 .278 .489 .000 .000 .115 .083 .148 .000 .000 11.88% 3.51% 20.65% .007 .014
Y12: Externalising .106 .001 .216 .046 .092 .031 .000 .062 .046 .092 9.02% 3.66% 15.04% .001 .002
Y12: Cognitive ability .029 -.004 .061 .095 .095 .029 -.004 .061 .095 .095 12.37% 6.80% 17.88% .000 .000
Y16: GCSE core subjects grade .027 .003 .051 .024 .073 .022 .002 .041 .024 .073 −0.89% −4.02% 2.54% .560 .560
Note. Est = Estimate, LB = Lower Bound 95% Confidence Interval, UB = Upper Bound 95% Confidence Interval. Significant (pBonferroni-Holm) effects are highlighted in green (increases) or red (decreases).
* Standardized mean difference = (Meanattrition + imputed - Meanoriginal) / SDoriginal
† Percentage change = (Varattrition + imputed - Varoriginal) / Varoriginal
‡ P values are Bonferroni-Holm adjusted within each attrition timepoint

Create Circular Heatmaps

Standardised Mean Differences

Code
# Create circular heatmap for SMD
heatmap_data_smd = imputed_comparisons_df %>%
  filter(stat == "smd") %>%
  mutate(
    outcome_labeled = rq1y_twin_labels_clean_extrashort[match(dataset, rq1y_twin_labels_clean_extrashort)],
    outcome_labeled = factor(outcome_labeled, levels = rq1y_twin_labels_clean_extrashort),
    Variables_wrapped = str_wrap(variable, width = 8),
    # Set fill value to NA if not significant, otherwise keep original sign
    fill_value = case_when(
      pval_adj < 0.05 ~ y,
      TRUE ~ NA
    )
  ) %>%
  # Use fixed order from rq6y_labels
  mutate(
    variable = factor(variable, levels = rq6y_labels),
    Variables_wrapped = factor(Variables_wrapped, levels = str_wrap(rq6y_labels, width = 8))
  )

 ggplot(heatmap_data_smd, aes(x = Variables_wrapped, y = outcome_labeled)) +
  geom_tile(aes(fill = fill_value),
            color = "black",
            size = 0.5) +
  geom_text(
    aes(label = ifelse(pval_adj < 0.05,
                       gsub("^(-?)0\\.", "\\1.", sprintf("%.3f", y)), ""),
        size = 4.5
    ),
    color = "black"
  ) +
  scale_size_identity() +
  coord_polar(theta = "x", start = 0, direction = 1, clip = "off") +
  scale_fill_gradient2(
    low = "#2166ac",
    mid = "white",
    high = "#d73027",
    midpoint = 0,
    limits = c(-0.3, 0.3),
    na.value = "white",
    guide = "none"
  ) +
  labs(
    subtitle = expression(frac(bar(X)[attrition + imputed] - bar(X)[original], SD[original])),
    tag = "A3"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 0, size = 13, color = "black"),
    axis.text.y = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    plot.title = element_blank(),
    plot.subtitle = element_text(hjust = 0.5, size = 17, margin = margin(b = -45)),
    plot.tag = element_text(hjust = 0, vjust = 0, size = 30, face = "bold"),
    plot.tag.position = "topleft"
  ) +
  {
    y_positions = seq_along(levels(heatmap_data_smd$outcome_labeled))
    labels = levels(heatmap_data_smd$outcome_labeled)

    lapply(seq_along(labels), function(i) {
      annotate("text", x = 0.5, y = y_positions[i],
               label = paste0(labels[i], " "), size = 5, hjust = 1, angle = 7, vjust = -.3, fontface = "bold")
    })
  }
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

Code
save_plot("8_rq6_circular_heatmap_smd", width = 8, height = 8, trim = TRUE)

📊 View Full Resolution SMD Heatmap (PNG)

Variance Differences

Code
# Create circular heatmap for variance differences
heatmap_data_var = imputed_comparisons_df %>%
  filter(stat == "var") %>%
  mutate(
    outcome_labeled = rq1y_twin_labels_clean_extrashort[match(dataset, rq1y_twin_labels_clean_extrashort)],
    outcome_labeled = factor(outcome_labeled, levels = rq1y_twin_labels_clean_extrashort),
    Variables_wrapped = str_wrap(variable, width = 8),
    # Set fill value to NA if not significant, otherwise keep original sign
    fill_value = case_when(
      pval_adj < 0.05 ~ y,
      TRUE ~ NA
    )
  ) %>%
  # Use fixed order from rq6y_labels
  mutate(
    variable = factor(variable, levels = rq6y_labels),
    Variables_wrapped = factor(Variables_wrapped, levels = str_wrap(rq6y_labels, width = 8))
  )

ggplot(heatmap_data_var, aes(x = Variables_wrapped, y = outcome_labeled)) +
  geom_tile(aes(fill = fill_value),
            color = "black",
            size = 0.5) +
  geom_text(
    aes(label = ifelse(pval_adj < 0.05,
                       paste0(sprintf("%.1f", y * 100), "%"), ""),
        size = 4.5
    ),
    color = "black"
  ) +
  scale_size_identity() +
  coord_polar(theta = "x", start = 0, direction = 1, clip = "off") +
  scale_fill_gradient2(
    low = "#2166ac",
    mid = "white",
    high = "#d73027",
    midpoint = 0,
    limits = c(-.126, .126),
    na.value = "white",
    guide = "none"
  ) +
  labs(
    subtitle = expression(frac(Var[attrition + imputed] - Var[original], Var[original]) ~ "\u00D7" ~ 100),
    tag = "B3"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 0, size = 13, color = "black"),
    axis.text.y = element_blank(),
    axis.title = element_blank(),
    panel.grid = element_blank(),
    panel.background = element_blank(),
    plot.background = element_blank(),
    plot.title = element_blank(),
    plot.subtitle = element_text(hjust = 0.5, size = 17, margin = margin(b = -45)),
    plot.tag = element_text(hjust = 0, vjust = 0, size = 30, face = "bold"),
    plot.tag.position = "topleft"
  ) +
  {
    y_positions = seq_along(levels(heatmap_data_var$outcome_labeled))
    labels = levels(heatmap_data_var$outcome_labeled)

    lapply(seq_along(labels), function(i) {
      annotate("text", x = 0.5, y = y_positions[i],
               label = paste0(labels[i], " "), size = 5, hjust = 1, angle = 7, vjust = -.3, fontface = "bold")
    })
  }

Code
save_plot("8_rq6_circular_heatmap_var", width = 8, height = 8, trim = TRUE)

📊 View Full Resolution Variance Heatmap (PNG)

Combine all pngs into one plot

Code 
Linking to ImageMagick 6.9.12.98
Enabled features: fontconfig, freetype, fftw, heic, lcms, pango, raw, webp, x11
Disabled features: cairo, ghostscript, rsvg
Using 16 threads
Code
gc()
           used  (Mb) gc trigger   (Mb)  max used  (Mb)
Ncells  5528539 295.3    8046225  429.8   8046225 429.8
Vcells 68297249 521.1  131988352 1007.0 124163771 947.3
Code
Sys.setenv(MAGICK_MEMORY_LIMIT = "8GiB")
Sys.setenv(MAGICK_MAP_LIMIT = "8GiB")
Sys.setenv(MAGICK_DISK_LIMIT = "16GiB")


# Read all circular plot images
img1 = image_read("plots/pngs/6_rq6_circular_heatmap_smd.png") %>% image_scale("33%")
img2 = image_read("plots/pngs/7_rq7_circular_heatmap_smd.png") %>% image_scale("33%")
img3 = image_read("plots/pngs/8_rq6_circular_heatmap_smd.png") %>% image_scale("33%")
img4 = image_read("plots/pngs/6_rq6_circular_heatmap_var.png") %>% image_scale("33%")
img5 = image_read("plots/pngs/7_rq7_circular_heatmap_var.png") %>% image_scale("33%")
img6 = image_read("plots/pngs/8_rq6_circular_heatmap_var.png") %>% image_scale("33%")


# Combine into rows
row1 = image_append(c(img1, img2, img3))  # SMD plots (A1, A2, A3)
row2 = image_append(c(img4, img5, img6))  # Variance plots (B1, B2, B3)

# Stack rows vertically
combined = image_append(c(row1, row2), stack = TRUE)

# Get dimensions for label positioning
info = image_info(combined)
row_height = info$height / 2

# Add row labels (adjust x and y coordinates as needed)
combined = combined %>%
  image_annotate(
    text = "Standardised Mean Difference",
    size = 50,
    degrees = -90,
    location = paste0("+50+", as.integer(row_height / 2)+450),
    color = "black",
    weight = 700
  ) %>%
  image_annotate(
    text = "Variance Percentage Change",
    size = 50,
    degrees = -90,
    location = paste0("+50+", as.integer(row_height * 1.5)+450),
    color = "black",
    weight = 700
  )

# Save combined image
image_write(combined, "plots/pngs/8_rq6_circular_plots_combined.png")
knitr::include_graphics("plots/pngs/8_rq6_circular_plots_combined.png")

📊 View plot

Correlations

Plot

Code
correlation_comparisons = list()

for(tp in seq_along(timepoint_names)){

  timepoint_name = timepoint_names[tp]

  # Get all imputations for this timepoint
  tp_pattern = paste0("^", timepoint_name)
  tp_indices = grep(tp_pattern, names(boot_compare_results))

  # Extract results for each imputation
  tp_results = boot_compare_results[tp_indices]

  # Extract correlation results across all imputations
  cor_data = lapply(tp_results, function(imp) {
    # imp$cor_resid is a list of bootstrap iterations
    sapply(imp$cor_resid, function(boot_iter) boot_iter)
  })

  # Combine all bootstrap samples across imputations
  all_cors = do.call(cbind, cor_data)

  # Calculate summary statistics
  cor_summary = apply(all_cors, 1, function(xx) .mean_qi_pd(xx))
  cor_summary = do.call(rbind.data.frame, cor_summary)
  cor_summary$x_var = x_var
  cor_summary$y_var = y_var
  cor_summary = cor_summary[cor_summary$x_var != cor_summary$y_var,]  # Remove diagonal
  cor_summary$pval_adj = stats::p.adjust(cor_summary$pval, method = "holm")

  correlation_comparisons[[tp]] = cor_summary
}

library(patchwork)

plots = lapply(1:length(correlation_comparisons), function(i) {
  plot_lower_triangular_matrix2(
    data = correlation_comparisons[[i]],
    variables = rq6y,
    labels = rq6y_labels,
    x_var_col = "x_var",
    y_var_col = "y_var",
    value_col = "y",
    p_col = "pval_adj",
    p_threshold = 0.05,
    show_values = TRUE,
    text_size = 3,
    title = rq1y_twin_labels_clean_extrashort[i],
    caption = NULL,
    method = "none"
  ) +
  labs(tag = paste0("C", i)) +
  theme(
    plot.title        = element_text(hjust=0),
    axis.text         = element_text(color = "black"),
    plot.background   = element_blank(),
    panel.background  = element_blank(),
    legend.background = element_blank(),
    plot.tag          = element_text(size = 14, face = "bold"),
    plot.tag.position = c(0, 1)
  )
})

patchwork::wrap_plots(
  plots,
  ncol = 3) +
  plot_annotation(
    title = expression(Cor[attritioned + imputed] ~ - ~ ~ Cor[original]),
    theme = theme(
      plot.title    = element_text(hjust = 0.5, size = 16, face = "bold"),
      plot.subtitle = element_text(hjust = 0.5, size = 13)
      )
  )

Code
save_plot(
  "8_rq6_correlations", width = 11, height = 3, trim = TRUE, trim_margin = 0)

📊 View Full Resolution Correlation Plots (PDF)

GT Results table

Code
cor_table = bind_rows(correlation_comparisons, .id = "timepoint") %>%
  mutate(
    dataset = rq1y_twin_labels_clean_extrashort[as.numeric(timepoint)],
    x_var_label = rq6y_labels[match(x_var, rq6y)],
    y_var_label = rq6y_labels[match(y_var, rq6y)],
    pair = paste(x_var_label, "×", y_var_label)
  ) %>%
  select(dataset, pair, y, ymin, ymax, pval, pval_adj) %>%
  gt(groupname_col = "dataset") %>%
  tab_options(
    table.width = px(800)
  ) %>%

  # Rename columns
  cols_label(
    pair = "Variable Pair",
    y = "Est",
    ymin = "LB",
    ymax = "UB",
    pval = md("p<sub>unadj</sub>"),
    pval_adj = md("p<sub>adj</sub>")
  ) %>%

  # Create column spanner
  tab_spanner(
    label = "Correlation Difference",
    columns = c(y, ymin, ymax, pval, pval_adj)
  ) %>%

  # Format numbers
  fmt(
    columns = c(y, ymin, ymax),
    fns = function(x) {gbtoolbox::apa_num(x, n_decimal_places = 3)}
  ) %>%
  fmt(
    columns = c(pval, pval_adj),
    fns = function(x) {gbtoolbox::apa_num(x, n_decimal_places = 3)}
  ) %>%

  # Styling - uniform font size
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_column_spanners()
  ) %>%
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_body()
  ) %>%
  tab_style(
    style = cell_text(size = px(10), align = "center"),
    locations = cells_row_groups()
  ) %>%
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_column_labels()
  ) %>%
  tab_style(
    style = cell_text(size = px(10)),
    locations = cells_footnotes()
  ) %>%

  # Highlight significant results - positive effects (light green)
  tab_style(
    style = cell_fill(color = "#d5e8d4"),
    locations = cells_body(
      columns = c(y, ymin, ymax, pval, pval_adj),
      rows = pval_adj < 0.05 & y > 0
    )
  ) %>%

  # Highlight significant results - negative effects (light red)
  tab_style(
    style = cell_fill(color = "#f8cecc"),
    locations = cells_body(
      columns = c(y, ymin, ymax, pval, pval_adj),
      rows = pval_adj < 0.05 & y < 0
    )
  ) %>%

  # Add footnotes
  tab_footnote(
    footnote = md("<em>Note.</em> Est = Estimate, LB = Lower Bound 95% Confidence Interval, UB = Upper Bound 95% Confidence Interval. Significant (p<sub>Bonferroni-Holm</sub>) effects are highlighted in green (increases) or red (decreases)."),
    placement = "right"
  ) %>%
  tab_footnote(
    footnote = "P values are Bonferroni-Holm adjusted within each attrition timepoint",
    locations = cells_column_labels(columns = pval_adj),
    placement = "right"
  ) %>%
  tab_footnote(
    footnote = md("Correlation difference = Cor<sub>attritioned + imputed</sub> - Cor<sub>original</sub>"),
    locations = cells_column_spanners(spanners = "Correlation Difference"),
    placement = "right"
  ) %>%
  opt_footnote_marks(marks = c("*", "†", "‡"))

cor_table
Variable Pair
Correlation Difference*
Est LB UB punadj padj†
Y21
Y12: Externalising × Y12: Depression (MFQ) .010 -.006 .028 .231 .886
Y12: Cognitive ability × Y12: Depression (MFQ) -.000 -.019 .019 .977 1.000
Y16: GCSE core subjects grade × Y12: Depression (MFQ) -.015 -.031 .001 .071 .428
Y12: Cognitive ability × Y12: Externalising .013 -.006 .032 .177 .886
Y16: GCSE core subjects grade × Y12: Externalising -.003 -.018 .013 .707 1.000
Y16: GCSE core subjects grade × Y12: Cognitive ability -.008 -.021 .004 .208 .886
Y26 (q'aire)
Y12: Externalising × Y12: Depression (MFQ) .003 -.017 .023 .762 1.000
Y12: Cognitive ability × Y12: Depression (MFQ) -.000 -.022 .022 .979 1.000
Y16: GCSE core subjects grade × Y12: Depression (MFQ) -.009 -.029 .011 .383 1.000
Y12: Cognitive ability × Y12: Externalising .009 -.013 .031 .453 1.000
Y16: GCSE core subjects grade × Y12: Externalising .007 -.011 .025 .462 1.000
Y16: GCSE core subjects grade × Y12: Cognitive ability -.017 -.032 -.001 .034 .206
Y26 (web test)
Y12: Externalising × Y12: Depression (MFQ) -.053 -.094 -.015 .007 .036
Y12: Cognitive ability × Y12: Depression (MFQ) .051 .006 .101 .025 .076
Y16: GCSE core subjects grade × Y12: Depression (MFQ) .013 -.019 .045 .428 .428
Y12: Cognitive ability × Y12: Externalising .045 -.001 .090 .054 .108
Y16: GCSE core subjects grade × Y12: Externalising .048 .013 .082 .009 .037
Y16: GCSE core subjects grade × Y12: Cognitive ability -.104 -.143 -.068 .000 .000
Note. Est = Estimate, LB = Lower Bound 95% Confidence Interval, UB = Upper Bound 95% Confidence Interval. Significant (pBonferroni-Holm) effects are highlighted in green (increases) or red (decreases).
* Correlation difference = Corattritioned + imputed - Cororiginal
† P values are Bonferroni-Holm adjusted within each attrition timepoint

Combine all correlation plots into one PNG

Code 
           used  (Mb) gc trigger (Mb)  max used  (Mb)
Ncells  5574031 297.7   10147732  542   8046225 429.8
Vcells 68501730 522.7  131988352 1007 124163771 947.3
Code
Sys.setenv(MAGICK_MEMORY_LIMIT = "8GiB")
Sys.setenv(MAGICK_MAP_LIMIT = "8GiB")
Sys.setenv(MAGICK_DISK_LIMIT = "16GiB")

# Read all correlation plot images
img1 = image_read("plots/pngs/6_rq6_correlations.png")
img2 = image_read("plots/pngs/7_rq7_correlations.png")
img3 = image_read("plots/pngs/8_rq6_correlations.png")

# Add vertical spacing (white border) between images
img1 = image_border(img1, "white", "0x20")  # Add 20px white space at bottom
img2 = image_border(img2, "white", "0x20")  # Add 20px white space at bottom

# Stack all plots vertically
combined = image_append(c(img1, img2, img3), stack = TRUE)

# Save combined image
image_write(combined, "plots/pngs/8_rq6_correlations_combined.png")
knitr::include_graphics("plots/pngs/8_rq6_correlations_combined.png")

ACE Comparison

Code
# Extract ACE estimates and differences from all datasets
ace_estimates_list = list()
ace_differences_list = list()

for(tp in seq_along(timepoint_names)){

  timepoint_name = timepoint_names[tp]

  # Get all imputations for this timepoint
  tp_pattern = paste0("^", timepoint_name)
  tp_indices = grep(tp_pattern, names(boot_compare_results))

  # Extract results for each imputation
  tp_results = boot_compare_results[tp_indices]

  # Process ACE results across all imputations
  ace_results = lapply(tp_results, function(imp) {
    lapply(imp$ace, function(x) x)
  })
  
  # Flatten and process
  ace_results = unlist(ace_results, recursive = FALSE)
  ace_results = lapply(ace_results, function(x) do.call(rbind, x))
  ace_results = do.call(rbind, ace_results) %>%
    filter(par %in% c("a","c", "e"))

  ace_results_summary = ace_results %>%
    group_by(par, name, sex, group) %>%
    summarise(
      .mean_qi_pd(value),
      .groups = "drop"
    )

  # Separate estimates and differences
  ace_estimates_list[[tp]] = ace_results_summary %>%
    filter(group != "diff")

  ace_differences_list[[tp]] = ace_results_summary %>%
    filter(group == "diff")

  # Add dataset identifier
  ace_estimates_list[[tp]]$dataset = rq1y_twin1[tp]
  ace_differences_list[[tp]]$dataset = rq1y_twin1[tp]
}
Warning: There was 1 warning in `summarise()`.
ℹ In argument: `.mean_qi_pd(value)`.
ℹ In group 39: `par = "c"`, `name = "lsdqext"`, `sex = "female"`, `group =
  "diff"`.
Caused by warning:
! pd-values smaller than 0.5 detected, indicating inconsistent direction
  of the probability mass. This usually happens when the parameters space
  is not continuous. Affected values are set to 1. See help('p_direction')
  for more info.
Code
# Combine into single dataframes
ace_estimates = do.call(rbind, ace_estimates_list)
ace_differences = do.call(rbind, ace_differences_list)

# Add readable dataset labels
ace_estimates$dataset_label = rq1y_twin_labels_clean[match(ace_estimates$dataset, rq1y_twin1)]
ace_differences$dataset_label = rq1y_twin_labels_clean[match(ace_differences$dataset, rq1y_twin1)]

Table - ACE estimates

Code
ace_estimates %>%
  filter(par %in% c("a", "c", "e")) %>%
  select(-starts_with("."), -n, -pd, -pval) %>%
  mutate(
    name_clean = rq6y_labels[match(name, rq6y_prefix)],
    group = case_when(
      group == "df1" ~ "Original",
      group == "df2" ~ "Imputed",
      TRUE ~ group
    )
  ) %>%
  pivot_wider(
    id_cols = c(dataset_label, name_clean, par, sex),
    names_from = group,
    values_from = c(y, ymin, ymax),
    names_sep = "_"
  ) %>%
  select(dataset_label, name_clean, par, sex,
         y_Original, ymin_Original, ymax_Original,
         y_Imputed, ymin_Imputed, ymax_Imputed) %>%
  arrange(sex, dataset_label, name_clean, par) %>%
  mutate(par = toupper(par)) %>%
  gt(groupname_col = "sex") %>%
  fmt_number(decimals = 3) %>%
  tab_row_group(
    label = "Male",
    rows = sex == "male"
  ) %>%
  tab_row_group(
    label = "Female",
    rows = sex == "female"
  ) %>%
  cols_hide(sex) %>%
  cols_label(
    dataset_label = "Dataset",
    name_clean = "Variable",
    par = "",
    y_Original = "Est",
    ymin_Original = "Lower",
    ymax_Original = "Upper",
    y_Imputed = "Est",
    ymin_Imputed = "Lower",
    ymax_Imputed = "Upper"
  ) %>%
  tab_spanner(
    label = "Original",
    columns = c(y_Original, ymin_Original, ymax_Original)
  ) %>%
  tab_spanner(
    label = "Attritioned + Imputed",
    columns = c(y_Imputed, ymin_Imputed, ymax_Imputed)
  ) %>%
  tab_options(
    table.font.size = px(9)
  )
Dataset Variable
Original
Attritioned + Imputed
Est Lower Upper Est Lower Upper
Female
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability A 0.473 0.336 0.614 0.428 0.312 0.549
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability C 0.209 0.079 0.332 0.271 0.159 0.379
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability E 0.318 0.282 0.357 0.301 0.271 0.332
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) A 0.211 0.000 0.443 0.219 0.001 0.435
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) C 0.236 0.042 0.410 0.205 0.021 0.377
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) E 0.553 0.477 0.628 0.576 0.508 0.643
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising A 0.513 0.422 0.568 0.493 0.396 0.549
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising C 0.009 0.000 0.083 0.012 0.000 0.091
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising E 0.478 0.430 0.529 0.494 0.450 0.541
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade A 0.582 0.493 0.676 0.544 0.461 0.631
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade C 0.307 0.215 0.392 0.347 0.263 0.426
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade E 0.112 0.098 0.126 0.109 0.096 0.122
Y26 (CATSLife web tests) Y12: Cognitive ability A 0.447 0.224 0.632 0.427 0.309 0.545
Y26 (CATSLife web tests) Y12: Cognitive ability C 0.152 0.000 0.342 0.272 0.162 0.381
Y26 (CATSLife web tests) Y12: Cognitive ability E 0.401 0.340 0.471 0.301 0.271 0.332
Y26 (CATSLife web tests) Y12: Depression (MFQ) A 0.235 0.000 0.494 0.220 0.004 0.431
Y26 (CATSLife web tests) Y12: Depression (MFQ) C 0.184 0.000 0.375 0.206 0.028 0.379
Y26 (CATSLife web tests) Y12: Depression (MFQ) E 0.582 0.476 0.684 0.574 0.505 0.643
Y26 (CATSLife web tests) Y12: Externalising A 0.355 0.151 0.490 0.494 0.394 0.549
Y26 (CATSLife web tests) Y12: Externalising C 0.067 0.000 0.232 0.013 0.000 0.093
Y26 (CATSLife web tests) Y12: Externalising E 0.578 0.503 0.657 0.494 0.449 0.541
Y26 (CATSLife web tests) Y16: GCSE core subjects grade A 0.703 0.580 0.848 0.544 0.461 0.631
Y26 (CATSLife web tests) Y16: GCSE core subjects grade C 0.182 0.042 0.301 0.347 0.262 0.427
Y26 (CATSLife web tests) Y16: GCSE core subjects grade E 0.115 0.098 0.135 0.109 0.096 0.123
Y26 (TEDS26 questionnaire) Y12: Cognitive ability A 0.451 0.308 0.597 0.426 0.310 0.547
Y26 (TEDS26 questionnaire) Y12: Cognitive ability C 0.224 0.089 0.352 0.272 0.161 0.379
Y26 (TEDS26 questionnaire) Y12: Cognitive ability E 0.325 0.287 0.368 0.301 0.272 0.332
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) A 0.239 0.000 0.475 0.216 0.000 0.431
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) C 0.207 0.007 0.397 0.208 0.026 0.381
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) E 0.554 0.472 0.633 0.576 0.509 0.643
Y26 (TEDS26 questionnaire) Y12: Externalising A 0.496 0.399 0.554 0.494 0.393 0.550
Y26 (TEDS26 questionnaire) Y12: Externalising C 0.008 0.000 0.083 0.012 0.000 0.092
Y26 (TEDS26 questionnaire) Y12: Externalising E 0.496 0.445 0.550 0.494 0.448 0.540
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade A 0.580 0.491 0.672 0.544 0.462 0.630
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade C 0.311 0.220 0.398 0.348 0.263 0.427
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade E 0.109 0.096 0.123 0.109 0.096 0.122
Male
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability A 0.472 0.274 0.674 0.369 0.240 0.502
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability C 0.243 0.062 0.410 0.374 0.251 0.490
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability E 0.285 0.234 0.347 0.257 0.223 0.295
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) A 0.358 0.134 0.487 0.396 0.260 0.470
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) C 0.049 0.000 0.224 0.014 0.000 0.115
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) E 0.593 0.507 0.684 0.591 0.529 0.653
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising A 0.447 0.268 0.554 0.545 0.431 0.605
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising C 0.043 0.000 0.190 0.019 0.000 0.113
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising E 0.509 0.442 0.584 0.436 0.393 0.483
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade A 0.630 0.499 0.773 0.615 0.516 0.717
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade C 0.211 0.075 0.334 0.258 0.157 0.352
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade E 0.159 0.136 0.185 0.127 0.111 0.145
Y26 (CATSLife web tests) Y12: Cognitive ability A 0.427 0.073 0.633 0.369 0.237 0.502
Y26 (CATSLife web tests) Y12: Cognitive ability C 0.035 0.000 0.277 0.374 0.248 0.493
Y26 (CATSLife web tests) Y12: Cognitive ability E 0.538 0.366 0.718 0.258 0.224 0.295
Y26 (CATSLife web tests) Y12: Depression (MFQ) A 0.155 0.000 0.367 0.395 0.264 0.469
Y26 (CATSLife web tests) Y12: Depression (MFQ) C 0.078 0.000 0.260 0.014 0.000 0.111
Y26 (CATSLife web tests) Y12: Depression (MFQ) E 0.767 0.615 0.898 0.591 0.530 0.654
Y26 (CATSLife web tests) Y12: Externalising A 0.271 0.000 0.466 0.545 0.430 0.605
Y26 (CATSLife web tests) Y12: Externalising C 0.053 0.000 0.257 0.019 0.000 0.115
Y26 (CATSLife web tests) Y12: Externalising E 0.676 0.532 0.813 0.436 0.393 0.482
Y26 (CATSLife web tests) Y16: GCSE core subjects grade A 0.729 0.566 0.806 0.616 0.520 0.716
Y26 (CATSLife web tests) Y16: GCSE core subjects grade C 0.017 0.000 0.162 0.256 0.158 0.349
Y26 (CATSLife web tests) Y16: GCSE core subjects grade E 0.254 0.194 0.327 0.127 0.111 0.145
Y26 (TEDS26 questionnaire) Y12: Cognitive ability A 0.477 0.241 0.689 0.370 0.239 0.503
Y26 (TEDS26 questionnaire) Y12: Cognitive ability C 0.178 0.000 0.375 0.373 0.248 0.490
Y26 (TEDS26 questionnaire) Y12: Cognitive ability E 0.345 0.272 0.427 0.257 0.223 0.294
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) A 0.236 0.000 0.436 0.396 0.264 0.470
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) C 0.120 0.000 0.321 0.014 0.000 0.113
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) E 0.643 0.550 0.734 0.590 0.529 0.653
Y26 (TEDS26 questionnaire) Y12: Externalising A 0.426 0.233 0.537 0.545 0.432 0.606
Y26 (TEDS26 questionnaire) Y12: Externalising C 0.041 0.000 0.199 0.019 0.000 0.112
Y26 (TEDS26 questionnaire) Y12: Externalising E 0.533 0.460 0.608 0.437 0.393 0.483
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade A 0.616 0.476 0.763 0.616 0.517 0.715
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade C 0.226 0.087 0.356 0.257 0.160 0.351
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade E 0.158 0.133 0.189 0.127 0.111 0.145

Table - Changes in ACE Estimates

Code
ace_differences_tbl = ace_differences  %>%
  group_by(dataset) %>%
  mutate(
    pval_adj = stats::p.adjust(pval, method = "holm")
  ) %>%
  ungroup()

ace_differences_tbl %>%
  select(-starts_with("."), -n, -pd) %>%
  mutate(name_clean = rq6y_labels[match(name, rq6y_prefix)]) %>%
  select(dataset_label, name_clean, par, sex, y, ymin, ymax, pval, pval_adj) %>%
  arrange(sex, dataset_label, name_clean, par) %>%
  mutate(par = toupper(par)) %>%
  gt(groupname_col = "sex") %>%
  fmt_number(decimals = 3) %>%
  tab_row_group(
    label = "Male",
    rows = sex == "male"
  ) %>%
  tab_row_group(
    label = "Female",
    rows = sex == "female"
  ) %>%
  cols_hide(sex) %>%
  cols_label(
    dataset_label = "Dataset",
    name_clean = "Variable",
    par = "",
    y = "Diff",
    ymin = "Lower",
    ymax = "Upper",
    pval = md("p<sub>unadj</sub>"),
    pval_adj = md("p<sub>adj</sub>")
  ) %>%
  tab_spanner(
    label = "95% CI",
    columns = c(ymin, ymax)
  ) %>%
  tab_style(
    style = cell_fill(color = "#d5e8d4"),
    locations = cells_body(
      columns = everything(),
      rows = pval_adj < 0.05 & y > 0
    )
  ) %>%
  tab_style(
    style = cell_fill(color = "#f8cecc"),
    locations = cells_body(
      columns = everything(),
      rows = pval_adj < 0.05 & y < 0
    )
  ) %>%
  tab_footnote(
    footnote = "Difference = (Attrition + Imputed) − Original. Adjusted P values are Bonferroni-Holm corrected within each sex, dataset, and parameter (A, C or E) group",
    locations = cells_column_labels(columns = pval_adj)
  ) %>%
  tab_options(
    table.font.size = px(9)
  )
Dataset Variable Diff
95% CI
punadj padj1
Lower Upper
Female
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability A 0.045 −0.061 0.154 0.408 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability C −0.062 −0.161 0.032 0.195 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability E 0.017 −0.009 0.047 0.234 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) A −0.008 −0.183 0.170 0.918 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) C 0.031 −0.114 0.172 0.662 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) E −0.023 −0.083 0.041 0.450 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising A 0.019 −0.047 0.099 0.486 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising C −0.004 −0.065 0.045 1.000 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising E −0.016 −0.056 0.025 0.440 1.000
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade A 0.038 −0.016 0.096 0.175 1.000
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade C −0.040 −0.097 0.010 0.121 1.000
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade E 0.003 −0.008 0.014 0.610 1.000
Y26 (CATSLife web tests) Y12: Cognitive ability A 0.020 −0.198 0.223 0.856 1.000
Y26 (CATSLife web tests) Y12: Cognitive ability C −0.119 −0.296 0.064 0.219 1.000
Y26 (CATSLife web tests) Y12: Cognitive ability E 0.100 0.039 0.169 0.000 0.005
Y26 (CATSLife web tests) Y12: Depression (MFQ) A 0.015 −0.281 0.306 0.922 1.000
Y26 (CATSLife web tests) Y12: Depression (MFQ) C −0.022 −0.246 0.212 0.845 1.000
Y26 (CATSLife web tests) Y12: Depression (MFQ) E 0.007 −0.102 0.116 0.896 1.000
Y26 (CATSLife web tests) Y12: Externalising A −0.138 −0.340 0.007 0.062 0.749
Y26 (CATSLife web tests) Y12: Externalising C 0.054 −0.050 0.216 0.499 1.000
Y26 (CATSLife web tests) Y12: Externalising E 0.084 0.007 0.164 0.029 0.403
Y26 (CATSLife web tests) Y16: GCSE core subjects grade A 0.159 0.037 0.299 0.008 0.140
Y26 (CATSLife web tests) Y16: GCSE core subjects grade C −0.165 −0.300 −0.049 0.004 0.080
Y26 (CATSLife web tests) Y16: GCSE core subjects grade E 0.006 −0.013 0.027 0.556 1.000
Y26 (TEDS26 questionnaire) Y12: Cognitive ability A 0.025 −0.093 0.143 0.685 1.000
Y26 (TEDS26 questionnaire) Y12: Cognitive ability C −0.049 −0.155 0.051 0.367 1.000
Y26 (TEDS26 questionnaire) Y12: Cognitive ability E 0.024 −0.008 0.060 0.154 1.000
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) A 0.023 −0.179 0.230 0.839 1.000
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) C 0.000 −0.163 0.159 0.999 1.000
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) E −0.022 −0.095 0.054 0.554 1.000
Y26 (TEDS26 questionnaire) Y12: Externalising A 0.002 −0.074 0.088 0.973 1.000
Y26 (TEDS26 questionnaire) Y12: Externalising C −0.004 −0.070 0.052 0.985 1.000
Y26 (TEDS26 questionnaire) Y12: Externalising E 0.002 −0.044 0.050 0.935 1.000
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade A 0.037 −0.027 0.105 0.273 1.000
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade C −0.037 −0.102 0.023 0.232 1.000
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade E 0.000 −0.012 0.013 0.966 1.000
Male
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability A 0.103 −0.086 0.297 0.283 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability C −0.131 −0.298 0.026 0.104 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Cognitive ability E 0.028 −0.025 0.089 0.322 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) A −0.038 −0.252 0.113 0.746 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) C 0.036 −0.071 0.204 0.721 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Depression (MFQ) E 0.002 −0.084 0.093 0.974 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising A −0.097 −0.264 0.030 0.123 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising C 0.024 −0.074 0.163 0.778 1.000
Y21 (TEDS21 phase-1 questionnaire) Y12: Externalising E 0.073 0.006 0.144 0.031 0.713
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade A 0.015 −0.101 0.140 0.820 1.000
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade C −0.047 −0.167 0.062 0.423 1.000
Y21 (TEDS21 phase-1 questionnaire) Y16: GCSE core subjects grade E 0.032 0.010 0.056 0.004 0.101
Y26 (CATSLife web tests) Y12: Cognitive ability A 0.058 −0.302 0.309 0.616 1.000
Y26 (CATSLife web tests) Y12: Cognitive ability C −0.338 −0.483 −0.090 0.011 0.187
Y26 (CATSLife web tests) Y12: Cognitive ability E 0.280 0.109 0.462 0.000 0.009
Y26 (CATSLife web tests) Y12: Depression (MFQ) A −0.240 −0.444 −0.006 0.044 0.572
Y26 (CATSLife web tests) Y12: Depression (MFQ) C 0.064 −0.082 0.254 0.620 1.000
Y26 (CATSLife web tests) Y12: Depression (MFQ) E 0.176 0.021 0.317 0.026 0.384
Y26 (CATSLife web tests) Y12: Externalising A −0.274 −0.565 −0.059 0.013 0.208
Y26 (CATSLife web tests) Y12: Externalising C 0.034 −0.096 0.244 0.939 1.000
Y26 (CATSLife web tests) Y12: Externalising E 0.240 0.092 0.379 0.000 0.009
Y26 (CATSLife web tests) Y16: GCSE core subjects grade A 0.113 −0.054 0.239 0.150 1.000
Y26 (CATSLife web tests) Y16: GCSE core subjects grade C −0.239 −0.343 −0.090 0.006 0.110
Y26 (CATSLife web tests) Y16: GCSE core subjects grade E 0.127 0.067 0.199 0.000 0.000
Y26 (TEDS26 questionnaire) Y12: Cognitive ability A 0.108 −0.136 0.332 0.367 1.000
Y26 (TEDS26 questionnaire) Y12: Cognitive ability C −0.195 −0.381 0.004 0.054 1.000
Y26 (TEDS26 questionnaire) Y12: Cognitive ability E 0.088 0.012 0.172 0.023 0.510
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) A −0.160 −0.403 0.045 0.176 1.000
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) C 0.107 −0.035 0.309 0.287 1.000
Y26 (TEDS26 questionnaire) Y12: Depression (MFQ) E 0.053 −0.041 0.148 0.294 1.000
Y26 (TEDS26 questionnaire) Y12: Externalising A −0.119 −0.310 0.020 0.089 1.000
Y26 (TEDS26 questionnaire) Y12: Externalising C 0.022 −0.082 0.175 0.897 1.000
Y26 (TEDS26 questionnaire) Y12: Externalising E 0.096 0.020 0.173 0.014 0.326
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade A 0.000 −0.129 0.134 0.986 1.000
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade C −0.030 −0.157 0.087 0.631 1.000
Y26 (TEDS26 questionnaire) Y16: GCSE core subjects grade E 0.030 0.005 0.061 0.015 0.354
1 Difference = (Attrition + Imputed) − Original. Adjusted P values are Bonferroni-Holm corrected within each sex, dataset, and parameter (A, C or E) group

ACE Estimates Plot

Code
ace_estimates %>%
  mutate(
    name_clean = rq6y_labels[match(name, rq6y_prefix)],
    name_clean = factor(name_clean, levels = rq6y_labels),
    group = case_when(
      group == "df1" ~ "O",
      group == "df2" ~ "I",
      TRUE ~ group
    ),
    group = factor(group, levels = c("I", "O")),
    sex = str_to_title(sex),
    par = toupper(par),
    dataset_label = factor(rq1y_twin_labels_clean_extrashort[match(dataset,rq1y_twin1)], levels = rq1y_twin_labels_clean_extrashort)
  ) %>%
  ggplot(aes(x = group, y = y, fill = par)) +
  geom_col(position = "stack", alpha = 1) +
  facet_grid(sex + dataset_label ~ name_clean, switch = "x") +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 0, hjust = .5, size = 8),
    strip.text = element_text(size = 8),
    strip.text.x = element_text(angle = 90, hjust = 1, vjust = .5),
    strip.text.y = element_text(angle = 0),
    strip.placement = "outside",
    panel.grid = element_blank(),
    plot.tag = element_text(hjust = 0, vjust = 0, size = 24, face = "bold"),
    plot.tag.position = "topleft"
  ) +
  labs(
    title = "Attritioned + Imputed vs Original",
    x = NULL,
    y = "Proportion",
    fill = "Component",
    tag = "C"
  ) +
  scale_y_continuous(breaks = c(.25,.50,.75, 1), labels = c(".25",".50",".75","1.00")) +
  scale_fill_manual(values = c("A" = "#d73027", "C" = "#fee08b", "E" = "#4575b4")) +
  coord_cartesian(ylim = c(0, 1))

Code
save_plot("8_rq6_ace_comparison", width = 5, height = 8)

📊 View Full Resolution ACE Comparison Plot (PNG)

Combine all pngs into one plot

Code 
           used  (Mb) gc trigger (Mb)  max used (Mb)
Ncells  6413296 342.6   10147732  542  10147732  542
Vcells 80794347 616.5  131988352 1007 131988339 1007
Code
Sys.setenv(MAGICK_MEMORY_LIMIT = "8GiB")
Sys.setenv(MAGICK_MAP_LIMIT = "8GiB")
Sys.setenv(MAGICK_DISK_LIMIT = "16GiB")


# Read all circular plot images
img1 = image_read("plots/pngs/6_rq6_ace_comparison.png") %>% image_scale("33%") %>% magick::image_crop("75%x100%+0+0")
img2 = image_read("plots/pngs/7_rq7_ace_comparison.png") %>% image_scale("33%") %>% magick::image_crop("75%x100%+0+0")
img3 = image_read("plots/pngs/8_rq6_ace_comparison.png") %>% image_scale("33%")


# Combine into rows
row1 = image_append(c(img1, img2, img3))  # SMD plots (A1, A2, A3)

# Save combined image
image_write(row1, "plots/pngs/8_rq6_ace_comparison_combined.png")
knitr::include_graphics("plots/pngs/8_rq6_ace_comparison_combined.png")

📊 View Full Resolution ACE Combined Comparison Plot (PNG)

Notes

  • Bootstrap results aggregated across 200 multiply imputed datasets
  • Multiple imputation accounts for uncertainty due to missing data
  • Confidence intervals reflect both bootstrap uncertainty and imputation uncertainty
  • P-values adjusted using Holm-Bonferroni method within each timepoint and metric combination